A  1-D  Granular  Gas  As  A  Knudsen  Gas1 

Patricio  Cordero*,  Jose  Miguel  Pasini*,  and  Rosa  Ramirez* 


*  Departamento  de  Fisica,  Universidad  de  Chile,  Santiago,  Chile 
1  CECAM,  ENS-Lyon,  46  Allee  d’ltalie,  669007,  France 


Abstract.  We  study  a  one-dimensional  system  of  inelastic  hard  particles  driven  by  thermal  walls.  For 
low  dissipation  in  collisions  the  system  may  be  treated  as  a  perturbation  on  a  non-interacting  system, 
being  thus  conceptually  similar  to  a  Knudsen  gas.  When  the  quasielastic  system  is  subjected  to  gravity 
a  kinetic  theory  approach  is  suitable,  giving  a  faithful  description  of  results  stemming  from  Newtonian 
molecular  dynamics  simulations.  Effects  derived  from  the  finite  size  of  the  particles  may  be  taken  into 
account  without  invalidating  the  rarefied  system  assumption.  When  the  system  is  not  subject  to  gravity 
and  held  between  two  thermal  walls,  the  boundary  conditions  determine  whether  the  system  will  be  fluidized 
or  clusterized.  The  transition  from  the  clustering  regime  to  the  fluidized  state  is  described. 


INTRODUCTION 

Granular  systems  are  many-particle  systems  characterized  by  dissipative  collisions.  When  subjected  to  a 
sufficiently  strong  excitation,  these  systems  may  behave  fluid- like  [1-4].  From  the  very  beginning  several 
authors  have  attempted  to  derive  hydrodynamic  equations  for  these  systems  [5,6].  A  more  recent  granular 
hydrodynamics  is  in  [7].  If  the  excitation  of  the  granular  system  is  through  a  permanent  energy  injection,  the 
fluid  may  stabilize  in  a  low-density  stationary  gaseous  state.  This  state  is  necessarily  a  nonequilibrium  state 
and  it  usually  is  inhomogeneous  as  well.  To  develop  the  basic  features  of  the  theory  of  gaseous  granular  systems 
we  restrict  the  analysis  to  the  simplifying  inelastic  hard-sphere  model  (IHS)  [5]. 

In  this  work  we  study  a  non-trivial  dissipative  system:  the  one  dimensional  system  of  inelastic  hard  spheres. 
This  simple  model  shares  some  of  the  phenomena  of  higher  dimensional  dissipative  systems,  such  as  transition 
from  a  condensed  to  a  fluid  state  [8]  and  inelastic  collapse  [9].  Hydrodynamics  and  kinetic  theory  have  been 
derived  for  such  systems,  when  free  of  clustering,  basically  from  the  associated  Boltzmann  equation  [10-12],  but 
there  is  a  crucial  fact  which  makes  the  1-D  case  special:  the  velocity  distribution  function  is  heavily  influenced 
by  the  boundary  conditions.  Therefore  some  of  the  1-D  results  we  present  here  could  qualitatively  be  extended 
to  higher  dimensions  for  Knudsen  gases,  or  even  standard  gases  near  the  boundaries. 

In  this  paper  we  study  a  1-D  system  of  N  point-like  particles  interacting  via  collisions  (with  and  without 
gravity  present)  that  conserve  momentum  but  dissipate  kinetic  energy.  We  do  theory  and  Newtonian  molecular 
dynamic  simulations.  The  inelastic  collision  rules  which  take  the  velocities  (ci,C2)  into  (c[ , c'2 )  are 

cj  =qa  +  (l-  q)c2 ,  c'2  =  qc2  +  (l-q)c1,  (1) 

where  0  <  q  <  1/2,  and  q  =  0  corresponds  to  the  perfectly  elastic  case.  The  usual  restitution  coefficient  is 
r  =  1  —  2q.  Our  study  was  originally  focused  on  describing  stationary  gaseous  states  of  the  system  but,  as  we 
will  see,  there  are  circumstances  when  the  system  develops  a  cluster  and  does  not  settle  to  a  stationary  state. 

If  particles  are  relabeled  after  each  collision,  so  that  in  (1),  c\  < — >  c'2  it  is  seen  that  the  limit  q  =  0  is 
tantamount  to  a  system  of  N  point  particles  that  pass  through  each  other  without  ever  interacting.  In  the 
case  q  ^  0  particles  pass  through  each  other  losing  some  energy  in  the  process,  as  if  moving  through  a  viscous 
background. 

The  limit  q  =  0  makes  the  system  trivial,  in  the  sense  that  particles  do  not  interact  with  each  other,  and 
Liouville’s  problem  becomes  a  one-free-particle  problem.  The  exact  kinetic  equation  becomes  Boltzmann’s 
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equation  itself  with  a  vanishing  collision  term  because  Liouville’s  distribution  is  the  simple  product  of  N  one- 
particle  distribution  functions.  The  inelastic  case,  with  q  sufficiently  small,  can  be  regarded  as  a  perturbation 
of  the  free  case.  It  has  been  seen  that  a  true  perturbation  theory  can  be  defined  when  the  parameter  is  not 
simply  q  but  the  product  qN  with  N  arbitrarily  large  [12] ,  all  this  assuming  that  no  clusters  are  formed. 

To  achieve  gaseous  states  a  permanent  energy  injection  is  needed,  which  we  have  chosen  to  be  a  thermal 
wall.  Namely,  a  particle  hitting  this  wall  bounces  back  as  if  coming  from  a  heat  bath  at  temperature  Twau. 

Our  work  has  focused  on  the  study  of  two  different  regimes.  On  one  hand,  stationary  states  can  be  found 
when  the  dissipation  rate  in  the  system  is  sufficiently  low.  In  this  case  we  show  that  a  theoretical  treatment 
based  on  Boltzmann’s  equation  is  suitable  for  describing  the  system.  On  the  other  hand,  at  high  dissipation 
rates,  clusters  may  form  in  the  system.  In  this  case,  the  kinetic  theory  approach  mentioned  is  not  suitable 
anymore,  and  simulational  studies  are  needed  for  developing  insight  into  this  regime.  As  we  will  see,  the 
dissipation  rate  depends  both  on  the  inelastic  coefficient  q  as  well  as  on  the  boundary  conditions  and  any 
external  field  applied  to  the  system. 

The  1-D  Boltzmann  equation  in  the  inelastic  point  particle  case  with  f  f(x,  c)  dx  dc  =  1  and  precollision 
velocities  cj  =  [(1  —  q)  c\  —  qc2\/(l  —  2 q)  and  c\  =  [(1  —  q)c2  —  qc\]/{l  —  2 q)  is 

+°tx  =  N  j  (  (1  -2q)2  ~  /(», *,«)/(*, <»,*))  (c-c2)dc2  (2) 

If  q  =  0  the  right  hand  side  can  be  shown  to  vanish.  Restricting  our  study  to  the  limit  AT  — ►  oo,  but  keeping 
qN  finite,  the  equation  reduces  to  [12] 

r\  r\  \  r\ 

gi+cfa.~9  f(x>  c> 0  =  <1N  -qc  [/0*h  c>  l)  M(x .  c>  *)]  (3) 

M(x,c,t)  =  J  \c  —  c'|(c  —  c')  f(x,c',t)dc' 

A  method  of  solving  this  equation  is  by  expanding  the  distribution  in  powers  of  qN  and  solving  consistently 
at  each  order. 

If  the  system  consists  of  rods  of  size  a  instead  of  point-like  particles  the  situation  changes  radically.  The 
systems  proceeds  via  inelastic  hard  collisions  (1)  exactly  through  the  same  set  of  velocity  states,  provided  the 
box  is  enlarged  from  L  to  LI  =  L  +  N  a.  In  this  respect  the  velocity  distribution  alone  is  the  same,  but  one 
cannot  describe  the  elastic  limit  as  equivalent  to  particles  passing  through  each  other  without  interacting.  If 
particles  are  still  relabeled  after  each  collision,  they  will  make  instantaneous  jumps  of  length  a  at  every  collision. 
If  again  the  passage  from  Liouville’s  equation  to  a  one  particle  equation  is  attempted,  the  approximate  equation 
will  have  a  non-trivial  collision  term.  Some  configurational  correlations  exist  and  they  have  to  be  described 
using  an  Enskog-type  equation  [14].  The  most  obvious  differences  are  obtained  looking  at  the  effects  of  a  on 
the  density  profile,  as  seen  later  on. 

Finally,  we  focus  our  study  of  the  transition  from  the  clustering  regime  to  fluidized  states  by  simulating  a 
system  of  point  particles  not  subject  to  gravity  between  two  thermal  walls  at  different  temperatures  Ijeft  and 
-fright  •  For  this  problem,  when 


A 


’-Light  Fjeft 

Fright  T  Figft 


(4) 


is  below  a  certain  threshold,  a  unique  cluster  is  formed  and  the  system  remains  in  a  time-dependent  state  that 
we  describe  further  on  and  of  which  a  detailed  account  is  given  elsewhere  [13]. 


SYSTEM  IN  AN  OPEN  BOX  UNDER  THE  ACTION  OF  GRAVITY 

Point-like  particles 

In  this  case  (grains  of  mass  m  =  1,  base  at  temperature  TbaSe  in  energy  units,  and  acceleration  of  gravity 
g)  there  is  no  cluster  formation  in  the  quasielastic  case,  and  Boltzmann’s  equation  (2)  can  be  applied  with 
no  restriction.  Equation  (3)  is  solved  perturbatively  and  it  gives  an  extraordinarily  faithful  description  of  the 
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FIGURE  1.  At  left  is  the  distribution  function  at  the  base  ( x  =  0)  normalized  to  1  evaluated  to  zeroth  (dot-dashed 
line),  first  (dashed  line)  and  second  order  (solid  line)  for  qN  =  0.2  and  N  =  200.  The  observed  distribution  differs 
considerably  from  the  zeroth  order  approximation.  In  spite  of  this,  the  second  order  approximation  already  matches 
extraordinarily  well  the  simulational  data  (open  circles).  The  solid  S  shaped  curve  underneath  represents  the  second 
order  correction.  Notice  that  this  correction  is  almost  zero  for  positive  velocities.  At  right  circles  show  in  a  logarithmic 
scale  the  density  observed  in  a  simulation  of  a  system  with  qN  =  0.1  and  N  =  200.  The  density  is  shown  evaluated  up 
to  zero  (dashed  line),  first  (dot-dash  line)  and  second  order  (solid  line).  It  can  be  seen  that  the  prediction  up  to  second 
order  reproduces  extraordinarily  well  the  density  at  heights  few  particles  reach. 


behavior  of  the  simulated  system  when  qN  is  not  too  large.  Figure  1  shows  the  velocity  distribution  function  at 
the  base.  Usually  predicting  the  shape  of  the  distribution  at  a.  wall  is  harder  than  getting  it  far  from  the  walls 
but  the  perturbation  method  is  powerful  enough.  The  abrupt  change  in  the  behavior  of  the  distribution  at  the 
base  about  c  =  0  can  be  understood  because  the  particles  near  the  base  and  going  up  are  dominated  by  the 
thermal  boundary  condition.  Hence  they  have  almost  no  information  from  the  bulk  and  they  are  described  by  a 
distribution  close  to  a  Maxwellian,  while  on  the  other  hand  the  particles  arriving  to  the  base  carry  the  influence 
of  their  passage  through  the  bulk.  These  two  distributions  meet  at  the  origin  of  phase-space  producing  the 
singular  behaviour  seen  at  left  in  Fig.  1. 

Defining  the  characteristic  height  hc  of  the  system  as  the  height  below  which,  on  the  average,  90%  of 
the  particles  are  found,  it  will  be  seen  that  the  present  theoretical  framework  gives  a  good  description  well 
beyond  hc.  This  height  depends  on  qN  but  in  the  range  0  <  qN  <  0.2  it  varies  from  2.7  to  1.65  (in  units  of 
Th&Se/ 9)  hence  taking  hc  2  as  the  characteristic  length  in  all  our  cases  is  reasonable.  Figure  1  also  shows,  at 
right,  the  density  of  the  system  up  to  a  height  x  =  6,  where  the  density  is  three  orders  of  magnitud  smaller 
than  at  the  base,  and  still  the  second  order  perturbative  solution  fits  well  the  observed  density. 

The  velocity  of  particles  falling  through  a  viscous  fluid  tends  to  a  limiting  value.  The  meaning  of  the 
maximum  that  the  velocity  distribution  has  for  c  <  0  at  x  =  0  seems  to  be  related  to  a  similar  phenomenon 
even  though  our  fluid  cannot  be  described  as  a  fluid  of  uniform  viscosity. 

We  used  a  system  of  N  =  200  particles  because,  as  suggested  at  left  in  Fig.  2,  which  shows  the  temperature 
profile  for  systems  of  different  size  (up  to  N  =  1000),  size  independence  is  already  reached  for  N  =  200  when 
qN  <  0.1.  Units  are  such  that  T  =  0.5  for  the  perfectly  elastic  system. 


A  gas  of  rods 

If  the  system  consists  of  rods  of  size  a  an  Enskog-type  kinetic  equation  has  to  be  used  to  consider  that  colliding 
particles  are  in  different  positions  and  that  there  are  excluded-volume  effects.  The  size  effects  (measured  with 
£  =  crN  g/Twaii)  compete  with  the  inelastic  effects  (measured  with  qN)  and  it  is  interesting  to  look  at  regimes 
where  they  are  somehow  comparable.  The  results  presented  in  Fig.  3  in  this  subsection  correspond  to  cases 
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FIGURE  2.  At  left  the  temperature  profile,  observed  in  systems  with  qN  =  0.01  (upper  set  of  points)  and  with 
qN  =  0.1  (lower  set  of  points)  for  systems  with  N  particles:  N  =  20  (rhombus),  N  =  200  (solid  circles),  and  N  =  1000 
(open  circles).  For  qN  =  0.1  the  profiles  with  N  =  200  and  N  =  1000  coincide  because  the  size  independent  behavior 
has  been  reached  while  the  points  for  N  =  20  show  a  higher  temperature.  For  qN  =  0.01  even  for  N  =  20  the  asymptotic 
behavior  of  the  temperature  has  been  reached  bacause  second  order  corrections  (where  the  size  appears  for  the  first 
time)  are  too  small.  At  right  the  permanent  energy  influx  through  the  base  is  dissipated  in  the  “volume” .  The  figure 
shows  the  heat  flux  Jq,  as  a  function  of  height  for  systems  with  N  =  200  and  qN  =  0.01  (squares),  qN  =  0.05  (triangles) 
and  qN  =  0.1  (circles).  The  solid  (dashed)  lines  give  the  respective  heat  fluxes  evaluated  up  to  second  (first)  order. 


( n(x)-n0(x) )  /  qN 


(  n(X)  -  n0(X) )  /  qN 


FIGURE  3.  At  left  the  profile  of  the  density  correction  -Aj  (n  —  no)  as  a  function  of  height  for  a  system  with  C  =  0.1 
and  qN  =  0.02.  The  dotted  line  considers  only  the  £  corrective  term  (up  to  second  order  in  £),  the  dashed  line  considers 
only  the  dissipative  correction  up  to  first  order  in  qN.  The  solid  line  shows  the  combined  effect  while  the  circles  are  the 
data  obtained  from  the  molecular  dynamic  simulations.  At  right  a  system  with  £  =  qN  =  0.02  with  both  corrections 
evaluated  up  to  first  order.  Notice  that  the  density  correction  cannot  be  too  large  near  the  base  because  of  the  size  of 
the  rods.  In  both  cases  theory  fits  excellently  the  simulational  data. 


in  which  the  effects  of  dissipation  were  evaluated  up  to  first  order  [14],  It  is  possible  to  check  that  the  pair 
correlation  function  has  a  non-trivial  structure  even  in  the  purely  elastic  case  [15]. 

Here  comes  into  play  one  distinctive  feature  of  the  one-dimensional  system:  taking  into  account  corrections 
due  to  the  finite  (but  small)  size  of  the  particles  is  not  inconsistent  with  the  rarefied  system  assumption. 
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FIGURE  4.  A  cluster  seen  in  phase  space.  The  system  consists  of  N  =  10000  particles,  qN  =  0.1  and  the  temperature 
at  both  walls  is  the  same.  This  particular  cluster  contains  approximately  half  the  particles  of  the  system  and  the  size 
of  the  cluster  is  about  2.5  xlO-5  the  size  of  the  linear  box. 


N  =  1000  qN  =  0.01  T^.T^.1 


FIGURE  5.  The  cluster  grows  continuously  with  time  until  it  hits  a  wall.  The  collision  provokes  the  evaporation  of 
part  of  its  particles,  this  vapor  pushes  the  cluster  away  from  the  wall  and  the  growing  process  is  resumed. 


Therefore  the  system  may  be  treated  as  a  perturbation  to  the  non-collisional  case.  By  expanding  the  collision 
term  in  Boltzmann’s  equation  in  a  double  series  in  qN  and  £  [14]  we  can  consider  to  a  different  extent  size  and 
dissipation  effects. 

Figure  3  compares  the  density  corrections  due  to  dissipation  and  to  the  finite  size  of  the  particles.  It  can 
be  seen  that  at  the  base  these  effects  compete:  dissipation  tends  to  increase  the  density,  while  the  size  of  the 
particles  tends  to  decrease  it.  The  combined  effect  matches  closely  the  results  obtained  from  simulations. 


THE  GAS  OF  INELASTIC  POINTS  IN  THE  ABSENCE  OF  GRAVITY 

Presence  of  a  cluster 

The  1-D  system  is  left  to  evolve  between  two  thermal  walls  at  temperatures  Teft  and  Tight-  We  choose 
Tieft  <  Tight  and  units  such  that  Teft  Tight  =  1-  The  parameter  A  defined  in  (4)  is  used  to  quantify  how  far 
from  the  symmetrical  case  (A  =  0)  is  the  system.  The  case  A  =  1  represents  an  infinitely  strong  temperature 
gradient. 

Our  simulations  show  that  for  A  =  0  a  cluster  unavoidably  forms  and  moves  at  random  about  the  center  of 
the  box.  If  A  0  but  small  a  cluster  still  forms,  but  it  roams  closer  to  the  colder  wall.  In  molecular  dynamics 
simulations  we  detect  clusters  by  looking  for  chains  of  particles  two  or  three  orders  of  magnitud  closer  than 
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jt  times  the  size  of  the  box.  Figure  4  shows  the  tiny  portion  of  phase  space  where  the  cluster  is  concentrated. 
The  spiral  shape  is  time-dependent,  as  the  cluster  continues  to  slowly  contract. 


N  =  1000  qN  =  0.01  T,.„  =  T,ight  =  1 


Time 

FIGURE  6.  At  top  the  cluster’s  position  as  a  function  of  time  sometime  before  hitting  a  wall  and  the  evolution 
afterwards.  At  bottom  the  evolution  of  the  size  of  the  cluster  in  the  same  time  lapse. 


N  =  1 000 


FIGURE  7.  Clustering  transition:  when  A  is  above  the  threshold  shown  in  the  figure  the  system  does  not  exhibit 
cluster  formation  and  is  completely  fluidized.  The  zone  between  the  two  curves  represents  a  regime  where  it  is  difficult 
to  decide  whether  clusters  will  form. 

Whenever  a  cluster  exists  it  keeps  growing  in  size.  This  makes  the  rest  of  the  system  more  rarefied  in  time, 
so  that  a  cluster-wall  collision  becomes  more  probable.  Figure  5  shows  the  size  of  the  cluster  (in  number  of 
particles)  as  a  function  of  time.  It  is  seen  that  there  are  sudden  size  reductions  which  correspond  to  the  instants 
when  the  cluster  hits  a  wall.  One  of  these  collisions  is  shown  in  detail  in  Fig.  6.  At  top  the  position  of  the 
cluster  is  seen  as  a  function  of  time  showing  the  exact  moment  when  the  collision  occurs.  At  bottom,  in  the 
same  figure,  it  is  seen  how  the  size  of  the  cluster  has  a  sudden  reduction  to  restart  the  growth  process.  Thus, 
not  only  is  the  system  highly  clumped,  but  it  is  in  a  non-steady  state  as  well. 
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N  =  1000  qN  =  0.1  Tm  =  0.5  T.ght  =  2 


N  =  1 000  qN  =  0.1  Tleft  =  0.84  T  =  1.190476 


FIGURE  8.  The  velocity  distribution  functions  at  the  left  wall  obtained  from  Newtonian  molecular  dynamics  (solid 
lines)  are  compared  with  the  distributions  obtained  from  DSMC  (dashed  lines)  for  a  system  of  N  =  1000  particles  and 
qN  =  0.1.  For  this  value  of  qN  a  cluster  forms  roughly  when  A  <  0.16.  At  left  the  case  A  =  0.60  which  is  far  from  the 
threshold  and  at  right  the  case  A  =  0.172. 


N  =  1000  qN  =  0.1  T}e1t  =0.5  T^.2 


N  =  1000  qN  =  0.1  T|afl  =  0.84  Tngh|  =  1.190476 


FIGURE  9.  Density  profiles  for  the  cases  depicted  in  Fig.  8. 


The  completely  fluidized  system 


By  increasing  A  the  frequency  of  cluster-wall  collisions  grows,  sometimes  obtaining  “windows”  in  which  the 
cluster  completely  evaporates.  By  further  increasing  A  these  windows  grow  until  a  point  is  reached  where  no 
cluster  is  formed.  In  this  fashion  a  totally  fluidized  stationary  state  is  achieved.  Figure  7  shows,  for  N  =  1000 
particles,  the  curve  that  separates  the  fluidized  from  the  clustering  regime.  When  A  is  below  but  close  to  the 
transition  curve  clustering  occurs  in  windows,  while  farther  from  the  transition  the  cluster  is  always  present. 

If  the  state  of  the  system  is  sufficiently  far  from  the  threshold  of  cluster  formation  it  may  be  described  with 
the  dissipative  Boltzmann  equation  (2).  But  since  the  velocity  distribution  obtained  from  MD  simulations 
exhibits  a  peculiar  non-Gaussian  feature  at  slow  velocities  which  resembles  the  singular  behavior  observed  when 
there  is  a  cluster  [13],  it  may  also  be  that  there  are  still  dominant  correlation  effects  that  make  Boltzmann’s 
equation  inapplicable.  To  discern  this  we  compare  results  from  Newtonian  molecular  dynamics  (MD)  with 
the  corresponding  results  arising  from  direct-simulation  Monte  Carlo  (DSMC),  which  neglects  correlations.  As 
seen  in  Figs.  8  and  9,  the  results  agree  very  well,  except  when  the  value  of  A  is  close  to  the  onset  of  clustering. 


881 


CONCLUSIONS 


We  have  established  that  a  quasielastic  1-D  system  under  gravity  may  be  thought  of  as  a  rarefied  gas:  the 
collision  term  of  the  associated  Boltzmann  equation  may  be  considered  small  and  the  boundary  conditions 
affect  the  state  of  the  system  throughout  the  volume.  Results  derived  using  perturbation  theory  match  well 
those  obtained  from  molecular  dynamics  simulations.  Effects  due  to  the  finite  size  of  the  particles  may  also  be 
considered  perturbatively,  consistently  with  the  Knudsen  gas  assumption. 

The  same  system  not  subject  to  gravity  was  also  studied.  Confined  between  two  walls  at  different  temper¬ 
atures,  fluidized  and  clustering  states  are  found.  The  distribution  function  exhibits  distinctive  non-Gaussian 
features.  When  the  system  is  kept  fluidized,  these  features  are  shown  not  to  arise  from  correlations  by  com¬ 
paring  MD  and  DSMC  results. 
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